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The conditions for the entry of vortices into type-II superconductors being in the Meissner and/or 
■ mixed state, are studied by both numerical and analytical solution of the Ginzburg-Landau equa- 

tions. A modulation instability of the superconducting order parameter is shown to occur when the 
kinematic momentum (or supervelocity) of the condensate reaches a threshold value at the super- 
conductor edge or surface. Due to this instability, vortices start to nucleate at the edge or surface 
, and then penetrate deeper into the sample. It is found that the presence of surface defects causes a 

' noticeable drop of the first penetration field and leads to a qualitative change of the magnetization 

curve. Based on these results a simple phenomenological model of the edge or surface barrier, taking 
^ , into account the effect of surface defects, is suggested. 
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I. INTRODUCTION 



Since the discovery of type-II superconductors until now, much interest was devoted to the effects that surface or 
edge barriers have on the magnetic and dissipative properties of these materials. The solution of the barrier problem is 
' closely connected with the condition for vortex entry into type-II superconductors. The problem of vortex penetration 
£h ' was widely considered in the literature [l|j|]. The majority of these studies [^-§]|J investigated the stability of the 



Meissner state by analyzing the linearized Ginzburg-Landau equations (GLE). Some studies to find the conditions 
for the entry of vortex semiloops into a superconductor B used the London model combined with the Gibbs free 
energy approach. A series of works was dedicated to direct numerical modeling of the vortex entry into type-II 



o 

superconductors ®^ by numerically solving the time-dependent Ginzburg-Landau equations . 

■ In the present paper a criterion of vortex entry into type-II superconductors, being in the Meissner and/or mixed 

state is formulated based on a numerical solution of the time-dependent Ginzburg-Landau equations (TDGLE), 

supplemented by an analytical study of the linearized GLE. We show that the obtained criterion is valid for both 

7-H ' bulk and thin-film type-II superconductors and allows, in particular, to estimate the first-vortex entry field H s for 

superconductors with an ideal surface. We also investigate the first-vortex entry field H en in the presence of defects 
"... 

i on the surface and its dependence on the size of these defects. It is found that H en may be much lower than H s but 
still exceeds the lower critical field H c \ (if we do not deal with granular superconductors). Besides, we investigate the 
effect of the surface defects on the magnetization curve of pin- free superconductors. On the basis of obtained results 
we suggest a phenomenological theory of edge or surface barriers for superconductors with both ideal and imperfect 
O ■ surface. 

The paper is organized as follows. In Sec. II the vortex entry conditions are studied in bulk superconductors on 
• i-H , the basis of a numerical (Sec. II A) and analytical (Sec. II B) solution of the GLE. In Sec. Ill the same problem is 
solved for thin superconducting films. In Sec. IV the effect of surface defects on the vortex entry and exit, and hence 
on the shape of the hysteretic magnetization curve, is considered for the example of bulk superconductors. On the 
basis of the obtained results a simple phenomenological model for an edge or surface barrier is suggested in Sec. V. 



II. BULK SUPERCONDUCTORS 

Let an infinite superconducting slab of width w (0 < y < w) be placed in a parallel magnetic field H — (0, 0, —H). 
The time-dependent GLE in this case are two-dimensional because the problem is homogeneous along the coordinate 
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z, 

^ = ~ [(-iV - A) 2 * + *(|*| 2 - 1)] + X , (1) 
BA 

— =Rc[r(-iVA)f]- K 2 VxVxA. (2) 

Here the length is scaled in units of the coherence length £(T), time in units of r = 47r er n A 2 (T)/c 2 (A being the London 
penetration length), the vector potential A in terms of $o/(27r£), where $o = ch/2e is the quantum of magnetic flux, 
<r„ is the normal-state conductivity, k — A/£ is the Ginzburg-Landau parameter, C is a relaxation constant [^0), x a 
random "force" which simulates fluctuations of the order parameter, [|| and Re means the real part. 

The boundary conditions for Eqs. (1,2) along the transverse coordinate y are: [V x A] z | a= o !to = —H and (—id/dy — 
A y )^\y—Q w — 0- F° r A and \& we assume periodic boundary conditions along the longitudinal coordinate x to model 
a sample of infinite length. 



A. Vortex entry condition - numerical calculations 

Equations (1,2) are solved numerically by employing a method similar to that described in [^J. Namely, we introduce 
a link variable Uj = exp(— i J Ajdxj) (xj = x,y) that allows to satisfy the gauge-invariance condition on the grid 
introduced by discretizing the variables x and y. To solve this system of time-dependent equations we use the Euler 
method. 

The parameters of the film are chosen as follows: width w = 25, length L = 50 (in units of £). Our numerical 
analysis shows that the results of our study are not sensitive with respect to an increase of L; the reason for this 
is discussed below in this Section. The Ginzburg-Landau parameter k = 2 and the constant C — 0.5 were chosen 
to minimize the computation time. From the microscopic theory it follows that C = 12 [ p~0|JTT| ] . however in JTl] ] 
arguments to prefer another selection of C are given. For each change of H in steps of height SH — 0.01 (in units 
of H C 2) we compute the solution of Eqs. (1,2). A series of iteration steps is performed until the magnetic moment 
becomes independent of time. The random "force" x 1S distributed uniformly within the interval — s < x < s; the 
obtained results did not depend on the noise amplitude s provided s <C 1. In our simulations we chose s = 0.001. 

During numerical solution of the Eqs. (1,2) we studied the condition at which the entry of vortices into the su- 
perconductor started. It turns out that the vortex nucleation is controlled by one single quantity: the supervelocity 
II = V</> — A (in dimcnsionlcss units) on the sample surface. As soon as this quantity at the sample edge, II° dg0 , 
reaches the critical value LI cr , a modulation instability of the superconducting order parameter occurs at the edge. 
The quantity LT cr ~ 1 practically does not depend on the external field for our selection of parameters. In what follows 
we consider vortex entry/exit conditions near the left (y = 0) edge of the sample; from the symmetry of the problem 
follows that these conditions near the opposite edge (y = w] are completely analogous. 

From our numerical solution of the GLE we calculate the magnetization curve of a superconductor (see Fig. 1). 
It is seen that the magnetization curve has a nonmonotonic serrated shape resulting from the threshold character of 
the vortex entry. Indeed, let at some value of the external magnetic field the inequality II° dgc > II cr be fulfilled on 
the edge of the superconductor. Then a chain of vortices with period I starts to nucleate at the slab edge, see Fig. 2. 
As shows our numerical calculation, the period of the vortex chain depends on the quantity All = |n x — LI cr |. The 
greater is |AII|, the less is Z, i.e. the more vortices per unit length enter a superconductor at one step. The entering 
vortices reduce the value of n° dgc , since the part of II X generated by an external magnetic field, 11^, and the part 
generated by the penetrated vortices, have opposite sign at the slab edge. The more vortices enter the slab at a 
single step of the magnetic field, the less becomes n° dgc . It is therefore necessary to further raise the magnetic field 
to fulfill the condition II° dgc > II cr again. Particularly large serrations of the magnetization can be obtained for a 
mesoscopic superconductor (of size comparable with the vortex core size £) , as it was discussed in BJl^] . As displays 
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our numerical analysis, an increase of the specimen size L only leads to small corrections to the obtained results; this 
is due to the fact that vortices enter the superconductor in form of a chain, the space period of which is much less 
than L, see Fig. 2. Besides, it turns out that the higher is An the less is the time of vortex entry. 

Table 1 gives the results of our analysis of the vortex entry into a zero-field-cooled superconductor. The quantity 
r is of the order of 10 -12 — 10~ 13 s and strongly depends on the parameters of the superconductor. 

Figures 3 and 4 show how II° dgo and j£ dgc depend on the external field H . It is easily seen that before reaching the 
value 0.41 (the entry field of the first vortices), the current density at the edge reaches the Ginzburg-Landau depairing 
current density jcL, however vortices do not enter the sample yet. With further increase of the applied field, the 
current density reaches a maximum value j > jcL and then slowly decreases. As compared to this j° d&c (H), the 
dependence n° dgo (_ff) has different shape. Namely, up to the barrier suppression field H s it increases almost linearly 
with H , but then it nearly saturates to some critical value H cr - Moreover, by comparing Figs. 1 and 3, one sees that 
the sharp change in n° dge (7J) results in a sharp change of the magnetization M. Besides, we see that the quantity 
n cr practically does not depend on the magnetic field H . Thus, Figs. 3 to 4 demonstrate the threshold character of 
vortex entry, and also the fact that the true criterion of vortex entry is the value of the supcrvelocity II° dgc at the 
edge, and not the current density j% d & e at the edge. 

B. Vortex entry criterion - analytical study 

To find the quantity H cr analytically, we consider the stationary Ginzburg-Landau equations 

AF/ k 2 + F(l - F 2 - IT 2 ) = , 



v x v x n + F 2 n = o, 

in which we introduced new variables F — \^\ and II. In this case it is more convenient to measure distances 
in units of A. The supervelocity |II| is in units of <I>o/27r£, and the magnetic field in units of &o/2n{;\ = ^/2H C , 
where H c = $o/ yf&TtXtl is the thermodynamic critical field Jll| . We seek the solution of these equations in the form 
F = Fq + f and II = Q + q, where / and q are small (|q| <C |Q|, / <C Fq) perturbations of the stationary values of 
the order parameter Fo and momentum Q. 

Disregarding terms of order higher than linear in q and / we obtain a system of three linear differential equations 
for the functions / and q, 

-4a/ + (1 - 3F 2 - Q 2 )/- 2F (Qq) = , 

Performing the Fourier transform with respect to the longitudinal coordinate x with wavevector k (physically this 
means we search for the solution as a chain of vortices with a period ~ we obtain the following system of 

equations: 

-TTT + K 1 ~ 3F o Q 2 fc/« 2 ) - 2F (Qq) = 0, (3) 

,2„ , , „ 1^2 



k 2 q y + ik^ ) + q y F 2 + 2fF a Q y = , (4) 
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— ^ + ik^- + q x Fg + 2fF Q x = 0, (5) 

supplemented by the boundary conditions: df j dy\ y —Q )W = 0, q y \ y =o,w = 0, and (ikq y — dq x /dy)\ y= o. w — 0. Hereafter, 
/ and q denote the Fourier-transforms of the perturbations / and q. 

As shows our numerical analysis of the equations (1,2) (with n = 2) the perturbation is localized at the slab surface 
over a length of about £. This feature will be analytically confirmed below for the case c> 1, at which perturbations 
are also localized near the surface over a length of several £ (i.e. much less than A). For such a distance from the 
surfaces one has Q y w as a result of our boundary condition Q° dgo = and the absence of vortices near to the 
surface. Besides, it is possible to expand Q x and Fq near the surfaces of the slab into a Taylor series (with respect to 
y) keeping only the linear terms. It is possible to take advantage of the relation Fq — 1 — Q 2 that follows from the 
GLE for the order parameter Fq ( neglecting the term with second derivative) . Inserting q y from Eq. (4) into Eq. (5) 
we obtain a system of two second-order differential equations, valid near the slab surface: 

\^T2 - 3F o - Qx~ ( k l^f) ~ 1F Q x q x = 0, (6) 



d 2 q x F 2 dq x d k 2 
dy 2 k 2 + F Q 2 dy dy k 2 + F 2 



q x F 2 + 2fF Q x = . (7) 



In the Meissner state Eqs. (6), (7) are valid not only on the slab surface, but also in the bulk, since in the Meissner 
state one has Q y = everywhere in the slab. Equations (6) and (7) have been derived and studied numerically in 
Refs. to check the stability of the Meissner state in bulk superconductors. 

We are interested in the solution of Eqs. (6) and (7) that vanishes deep inside the slab, choosing the surface at 
y = 0. In what follows we shall concentrate on the physically most interesting limit k> 1. This allows us to neglect 
the term with the first derivative in Eq. (7) since the function k 2 / (k 2 + Fq) depends on y rather smoothly (in the limit 
fc«K>l and Fq < 1) . In addition, the choice k> 1 allows us to disregard the term with the second derivative in 
Eq. (6) at the critical wavenumber k = k c (see below). As a result we arrive at the Airy equation for the perturbation 

q x , 

^ + <&(& + ay) = 0. (8) 
dy 2 

It is worth noting here that Eq. (8) with proper boundary conditions is applicable to study the stability of both 
Meissner state and mixed state in bulk superconductors with respect to the vortex entry . In the latter case the 
validity region of Eq. (8) is yo <C w v f, where yo is the size of perturbation in the y-direction and w v f is the width of 
the vortex-free zone near the superconductor surface. In Eq. (8) the coefficients a and b are expressed through the 
constants P and R, which appear while expanding the order parameter Fq sw 1 — P and the vector potential related 
quantity Q 2 as P + Ry near the superconductor edges, 

b = -(1 + k 2 - P)(GP - 2 - {k/n) 2 )/{2P - 2 - {k/n) 2 ) , 



a = -4R(1 + k 2 - P)[3P - 2 - (k/n) 2 )/(2P - 2 - {k/n) 2 ) 2 . 

The coefficient R = — 2v r PH c< ig C is expressed through the magnetic field near the edge, -Hedge, which for bulk 
superconductors is equal to the external field H; the absence of linear terms in the expression for F^(y) follows from 
the boundary condition dFo/dy\ y= o = 0. The solution of the Eq. (8) is the Airy function 

q x = D^[-a 1 ^(y + b/a)}, 

where D is an integration constant. 
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From the boundary condition dq x /dy\ y= o = (derived from the conventional condition q y \ y =o — 0) follows the 
simple relation 

b = - 7 ia 2/3 , (9) 

where 71 w — f .02 is the first zero of the first derivative of the Airy function. Equation (9) defines the value of the 
supervelocity II on the slab surface for which a periodic, spatially modulated disturbance q x exists for a given 
value of the perturbation wave vector k. Inserting the expressions for a and b into Eq. (9) we obtain the dispersion 
relation in the form 

(P - 1 - k 2 ) (6P - 2 - {k/nff (2P-2- (k/n) 2 ) = gP (3P - 2 - {k/n) 2 f (10) 
that contains only one field-dependent control parameter 

3 = 64| 7l | 3 J ff c 2 dgc . 

The given equation relates the value II on the slab surface to the wavenumber k = k c for which a nontrivial solution 
q x becomes possible. 

Numerical analysis of Eq. (10) shows that the minimum (but non-zero) value of P(fc) is reached for k = k c . In 
Fig. 5 the dependence H x dgc (k) = \J P(k) is presented for different values of k and H. It is important that the 
value II cr = \J P(k — k c ) for k = 2 is practically independent of the magnetic field H. Indeed, for H = QA1H C 2 
[g(H = 0AH c2 ) w 45] H cr is equal to ~ 0.78, while for H ~ H c2 [g{H ~ H c2 ) « 280] H cr ~ 0.79 differs only slightly. 
Numerical results confirming this independence of H cr from 7J are presented in Fig. 3. A slight quantitative difference 
is possibly due to the fact that Eq. (10) was derived in the limit k ^> 1. Note also that for the latter case the difference 
between the critical values II cr at the field of first vortex entry (i.e. at H = H s w H c and at the field near H C 2 
becomes more appreciable (0.66 and 0.81, respectively, for k = 100); however the difference in AII cr does not exceed 
20%. 

Our analysis of the stability of the Meissner state for the case k = 2, performed on the basis of the linear Eqs. (6) 
and (7), gives a critical value II cr = 1.02 (for H = 0.41if C 2), which is quite close to the value H cr = 0.97 produced 
by direct solution of the nonlinear GL equations (1) and (2), but differs appreciably from the result of Eq. (10): 
n cr = 0.78. A similar comparison performed for the case k = 5 gives a more satisfactory coincidence of the II cr values 
(0.8, 0.85 and 0.75, respectively), which were obtained by solving Eqs. (l)-(2), (6)-(7) and (10), respectively. 

It is interesting to note that for the case k 3> 1, the perturbation is localized at the superconductor surface; this is 
reflected by the exponential decay of the Airy function <i>(0) in the region 9 > 0. The estimate for the characteristic 
scale yo of the perturbation decay deep into the film yields yo ~ b~ x / 2 . The numerical solution of (10) shows that for 
large values of k one has b(k c ) ss k 2 « 0.1k 2 (see Fig. 5); therefore, we obtain finally that the perturbation is localized 
on a scale of several £ along the y-direction. Besides, for k ^> 1, the important assumption d 2 f/dy 2 <C fn 2 is justified 
for k = k c w (0.3 — 0.4)k, which permits to simplify Eq. (6). 

C. Vortex exit 

Let us discuss now the conditions for the vortex exit which takes place when the external field is decreased. As 
can be seen from Figs. 1 and 3, the dependences M(H) and II° dge (if) are nearly linear functions of the (initially) 
decreasing magnetic field. Physically this is explained by the fact that during the initial stage of the field decrease the 
vortices do not exit from the superconductor. At some value of the magnetic field (in our case at H = 0.367J C 2) the 
vortices start to exit from the sample. At H < 0.36iJ C 2 both dependences M(H) and n° dgc (i?) become nonmonotonic 
because not one but several vortices leave the sample at once. Our study shows that the process of the vortex exit 
is controlled by the behavior of II° dgc : as soon as n® dge changes sign with the decrease of H, the vortices leave the 
superconductor. In the opposite case the vortices do not exit. 
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At first sight Fig. 3 seems not to confirm the above statement; this is explained by the fact that only stationary 
values of H x dgc are shown there. It turns out that during the initial decrease of the magnetic field n° dgG also decreases. 
For example, at H = 0.37H C 2 one has n° dg0 = 0.04; when the magnetic field decreases down to H — Q.36H C 2, n° dgc 
first becomes negative, and then, after exit of several vortices, reaches the stationary value 0.17. This process repeats 
at lower values of the magnetic field, e.g. at H/H C 2 = 0.25,0.16,0.09 for given sample/material parameters. 

We should emphasize that when the field decreases below some value H = H* (in our case H* = 0A1H C 2), the 
magnetization and magnetic field coincide in sign. Such a "paramagnetic" response, resembling the result obtained 
by Jl4| for the bulk case, is essentially due to surface pinning, which is a generic feature of type-II superconductors 
with an ideal surface. 



III. THIN FILMS 

Extending the method developed for bulk superconductors, we shall find below the criterion for vortex entry into a 
thin (d < A) superconducting film placed in a perpendicular magnetic field. In this geometry II is practically indepen- 
dent of z; together with the boundary conditions d\'$'\/dz\ z= ± c i/2 — this gives a z-independent order parameter 
By averaging equation (1) over the coordinate z one arrives at the two-dimensional GL equation. Using the method 
suggested by [[l5| for thin films, we write the GLE for II valid in all space as 

V x V x n = -F 2 Tl5(z)9(y)d(w - y) . (11) 

The equation to find the order parameter, valid for superconducting domain z = 0, < y < w, — oo < x < +oo, reads 

K~gAF + F(l - F 2 - II 2 ) = , (12) 

where K c ff = A G ff/£, A G fr = A 2 /d, and the distance now is measured in units of A G fr, II in units of $ /27r£, and the 
magnetic field in units of <&o/27r£A e ff ■ Similarly to the case of bulk superconductors, we linearize Eqs. (11) and (12) 
with respect to small (|q| -C |Q|, / -C Fq) perturbations of the stationary values of the order parameter F and 
momentum Q. As a result we obtain the combined equations for / = F — Fq and q = II — Q, 

K~ 2 Af + f(l - 3F 2 - Q 2 ) - 2F Qq = , (13) 

V x V x q = -(F 2 q + 2fF Q)S(z)9(y)9(w - y) , (14) 

We now perform the Fourier transform with respect to the coordinates x,z and account for the condition q z (z = 
0) = that follows from the condition j z = in the superconductor. It is easy to show, that the magnetic field 
generated by the perturbation h = V x q can be neglected (this applies at k — k c as displays further analysis). As a 
result we obtain the following combined equations for the Fourier transforms q x (k,y, z = 0) and f(k,y, z = 0) 

-^ r ^ 2 - + f(l-3F 2 -Q 2 (k/ KcS f)-2F {Q x q x + Q y q y ) = 0, (15) 

d 2 q x F n 2 dq x d F n 2 + k 2 d , , , . . 

- W - tTy^W- + d-y^ F ^ + ^ + 2 ^ - ° • ( 16 > 

In full analogy with the case of bulk superconductors, the "most dangerous" perturbations of the mixed state are 
localized near the film edges, provided K e ff 3> 1- Therefore, following our above analysis for bulk samples, we neglect 
(in the limit re e ff ^ 1) the component Q y at the edge as well as the term with the second derivative in Eq. (15). It 
is also useful to expand Fq and Q x in a Taylor series in y near the film edge. Substituting then / from Eq. (15) into 
Eq. (16) we obtain the Airy equation for q x valid near the film edges, 

d 2 q x 



dy 2 



+ q x (b f + a f y) = 0. (17) 
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The Airy-type equation was derived for the first time by [H when authors analyzed the stability of the Meissner state 
of wide (w ^> X 2 /d) thin strips carrying a current. A similar equation was later analyzed by |J when investigating the 
stability of the Meissner state for narrow (w <C X 2 /d) thin strips. We should emphasize that our approach extends 
the above studies to the more general problem of the stability of both the Meissner state and the mixed state with 
respect to vortex entry. 

In Eq. (17) the coefficients a,f and bf are expressed through the stationary values P and R, which are the coefficients 
of the Taylor expansion of Q 2 and Fq near the edges: Q 2 » P + Ry and F$ « 1 — P (R = —2VPH e dge), 

b f = -k 2 (6P - 2 - (k/ KcS ) 2 ) /(2P-2- (k/n cS ) 2 ) , 

a f = ~ARk 2 (3P - 2 - {k/n cS ) 2 ) /{2P-2- (k/ KcS ) 2 ) 2 . 

The expressions for the coefficients a/ and bf are equivalent (if we rescale y in units of A) to the corresponding 
expressions for bulk superconductors (as was already mentioned, Eqs. (8) and (17) are valid only for k» 1, Keg 3> 1, 
and the relevant values of the critical wavenumber k c satisfy the inequality k c 3> 1 in both cases). Therefore, all 
further analysis is completely analogous to the above analysis, with mere replacement of H by -Hedge in the coefficient 
R. 

Thus, for both superconducting films and bulk superconductors, the criterion for vortex entry is met when the 
kinematic momentum II at the sample edge or surface reaches the critical value II cr defined by Eq. (10) (in the limit 
k, n e g 3> 1). The quantity II cr in general only slightly depends on the external field H and GL parameter k. The 
minimum value of IT cr = 1/Vo is reached at g = (which corresponds to the GL depairing current density jgl), and 
the maximum value H cr = \/2/3 at g 1. It is worth noting that the critical value of II obtained in ^] for the first 
chain of entering vortices is close to our value. 



IV. EFFECT OF SURFACE DEFECTS ON THE CONDITION OF VORTEX ENTRY AND EXIT 

In the above study we concentrated on the vortex entry and exit conditions for superconductors with an ideal surface 
or edge. However, it is evident that surface defects, being an inevitable property of real superconductors, should 
strongly affect their magnetic properties. Therefore we study the effect of surface defects on the vortex entry and exit, 
and thus on the magnetic characteristics of superconductors. Our investigation is performed for bulk superconductors, 
and it is additionally supposed that the defects are homogeneously distributed along the coordinate z. On one side, 
this mimics columnar-like defects, on the other side it allows us to treat the problem as two dimensional. 

We consider defects of two basic types. The type-I defect represents an inclusion on the superconductor surface 
of a domain of a phase possessing a critical temperature T* different from the bulk T c . This defect is simulated by 
introducing into the GL equation for the order parameter \t a function p(r) ~ (T* — T)/ (T c — T) M that characterizes 
the defect strength, 

^ = -^[(-iV - A) 2 * + *(|*| 2 - p(r))j + x . (18) 

Note that p can vary from — oo up to +oo (inside the defect area), depending on the relative value of T,T C , and T*; 
in the defect-free area one has p(r) = 1. 

The type-II defect simulates the roughness of the surface (geometrical defects). On each face of the geometrical 
defect the following boundary conditions are used, V x A| s = —H and (— iV — A)W| n = 0. For both defect types we 
shall concentrate on the analysis of square-shaped defects. Besides this, we suppose that each slab surface possesses 
one identical defect per length L along the x direction, see Fig. 6. To be specific, in what follows we consider mainly 
the case k = 2. 
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A. Inclusion defects 



Let us consider first the influence of typc-f (inclusion) defects on the vortex entry and exit. Let p(r) = in the area 
of a defect, and the superconductor parameters are the same as for a superconductor with an ideal surface. Later on 
we shall discuss other functions p(r). 

Our computation shows that the value of the supervelocity II in the close vicinity of the defect corners considerably 
exceeds its value on the slab boundary far from the defect. Unfortunately the linearized analysis of GLE presented in 
section II is not applicable near the defect. Nevertheless, it is reasonable to expect that there exists a threshold (critical) 
value 11^, typical for a specific defect, at which the vortex nucleates and penetrates deep into superconductor. Our 
numerical solution of the TDGLE reveals that 11^ in fact only slightly depends on the physical properties and shape 
of the defect . Indeed, the vortices nucleate at the places where the local value of II reaches its maximum value (which 
is close to one in dimcnsionlcss units). In the case of a larger defect this maximum occurs on the corners of the defect. 
Note that the local current density in those places is significantly less (by a factor 2 or 3) than the depairing current, 
due to the local suppression of the order parameter. Thus, for a defective surface, the vortices start to move deep 
into the sample at the "real entry field" H en , which is lower than the vortex entry field H s for a superconductor with 
an ideal surface. However, as can be clearly seen from table II, for defects of size exceeding certain saturation size 
3£ x 3£, the field H en ceases to decrease. 

Moreover, our study shows that when the size of a defect exceeds 3£, the defect starts to pin "vortices" already 
penetrated inside the defect. For example, taking a defect of size 4£ x 4£ one finds that the critical value of the local 
supervelocity H cr at the point of fluxoid entry (the point A in Fig. 6) equals 0.94, while the local current density at 
vortex entry is only 0.17j s . The analogous data for defects of size 6£ x 6£ are as follows: 1.06 for II cr and 0.08j s 
for j at the same point A in Fig. 6. It is worth noting here that the supercurrent density is still much less than the 
depairing limit, since the order parameter is strongly suppressed inside the defect. 

Our above results clearly show that while the supervelocity at the moment of vortex entry has a rather universal 
value, the vortex-entry current density at that moment is extremely sensitive to the defect size and shape. This 
finding demonstrates the nonuniversality of the current density at vortex entry. We emphasize that the vortices which 
have penetrated into the defect are rather of Josephson nature and hereafter will be referred to as "J- vortices" . Such 
J-vortices are characterized mainly by the phase variable since the modulus of the condensate wave-function ^ is 
strongly suppressed (\^\ <C 1) inside the defect. The larger is the defect, the larger is the number of J-vortices that 
can penetrate into it. Figure 7 shows the longitudinal average of the supervelocity II X at the slab surface as a function 
of the applied magnetic field for defects of different sizes. One can see that during further increase of the field the 
J-vortices enter the defect one-by-one. Typical jumps Ail° dgc of the quantity II° dgc are visible; each jump corresponds 
to the entry of one J-vortex into the defect. Note that the height of a jump practically does not change with the 
increase of the magnetic field. For greater defect size more jumps occur when the magnetic field increases from zero 
up to the field H en for vortex entry deep into the superconductor. 

The presence of surface defects, besides of lowering the field of vortex entry deep into the superconductor, has also 
a considerable influence on the shape of the magnetization curve. In Fig. 1 (see curve 2a) the magnetization curve 
of a superconductor with defects of size 2£ x 2£ is presented. It is seen that the presence of type-I defects leads to 
smoothing of the dependence M{H) as well as of the dependences il^. dge (iJ) and j| dge (i?) shown in Figs. 3 to 4. 
Indeed, at the selected defect configuration (one defect per length L on each edge) only two vortices can enter a slab 
at a single step of the field increase. The penetrated vortices almost immediately reduce IT near the defect below the 
critical value. As a result, the magnetization curve becomes smooth. 

On the decreasing branch of the applied magnetic field, from its maximum value down to zero, the dependences 
M(H), II® dge (iJ), and j% dse {H) in the presence of surface defects qualitatively differ from the respective dependences 
in the case of an ideal surface. This is explained by the fact that the vortices start to exit a slab practically at 
once when the decrease of the magnetic field starts. Note that both entry and exit occur through surface defects. 
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Therefore, in the presence of surface defects M(H) becomes smoother, since a rather small number of vortices exit 
at each stepwise decrease of the field. However, the descending branch of the curve M(H) (field decreases) does not 
coincide with the ascending branch of the curve M (H) (field increases) . This shows that in the presence of surface 
defects the curve M(H) exhibits a hysteresis. Such an irreversible behavior of superconductors with a perfect and 
imperfect surface is in qualitatively good agreement with the results of Ref. | fL6| , where superconductors with both 
smooth and artificially rough surface were studied. 

Except for the above differences, there is also an important feature that is generic for both superconductors with an 
ideal surface and with an imperfect surface. Namely, we find that with increasing magnetic field the quantity n° dgc 
practically does not depend on the magnetic field in the field range H > H s (for superconductors with ideal surface) 
or in the range H > H en (for superconductors with surface defects). 

We wish to emphasize that for a slab with a defect of size 6£ x 6£ the field H en exceeds the lower critical field H c i . 
Indeed, the H c \ value calculated by numerical solution of the GLE for k = 2 is 0.19iJ C 2, while the vortex entry field 
H en from our result is 0.28iJ C 2. The test calculation with k — 5 for a slab of width w = 50£ = 10A and defect size 
6£ x 6£ has shown that in this case the vortex moves deep into the slab at a field appreciably exceeding H c \. Indeed 
for k — 5 the vortex entry field of a superconductor with defect-free surfaces, H s is 0.145iJ C 2; with surface defects of 
the given size H en is 0.085ff C 2, while one has H c \ — 0.044i? C 2. 

As already mentioned, all these results correspond to defects with p(v) = 0. It turns out that qualitatively similar 
results hold for arbitrary superconducting defect with p(r) > as well. For example, for a defect of size 6£ x 6£ with 
the choice p(r) = 0.5 one finds that H en equals 0.3H C 2, while for p(r) — one has H en — 0.28i? C 2. Note that on 
reaching this field the defect contains two vortices, while for p(r) = there are three vortices inside the defect. As 
regards defects with negative p(r) < (nonsuperconducting inclusion) we observe a tendency for saturation. Taking, 
for example, p(r) = —2.0 one finds that H en — Q.3H C 2 with two vortices being pinned by the defect. It turns out that 
with further decrease of p(r) the vortex entry field does not vary and two vortices remain inside the defect. 

B. Geometrical defects versus inclusion defects 

Finally, we shall consider surface roughness defects (of type-II in our classification) having the same shape and 
sizes as considered above. Naturally, vortices cannot enter such a defect since by definition in this place the order 
parameter is zero because of absence of free electrons. This situation permits to study only the field H en for vortex 
entry deep into the superconductor. In table III the dependence of on the size of the type-II defect is given. 

Comparing this with table II we see that for defects of equal sizes the entry field H* n is lower than it was for defects 
of type-II, -ffen, i- e -' ^'n < Hen- This means that inclusions of another phase depress H en much stronger than does 
surface roughness. For example, for a defect of size 6£ x 6£, equals 0.34i? c2 , while for a defect of the first type with 
p(r) — the entry field H.\ n equals 0.28_ff C 2. The following arguments may be formulated to explain the dominance 
of the inclusion defects in depressing H en . First, for the case of type-I defects (inclusions of a material with lower T c : 
p(r) < 1) up to the field H en there are at least two vortices inside sufficiently large defects (of size exceeding 4£ x 4£). 
These vortices, being localized near the corners of the defect, provide an additional (in fact, positive) contribution 
to the value of II inside the corners. This results in a significant decrease of H en , compare tables II and III above. 
Indeed, in such a situation a lower external field is needed to accelerate the superconducting condensate up to the 
critical momentum ET*:^ that has to be reached at the corners of the defect. 

A second reason for the higher sensitivity of H en with respect to inclusion defects is related to the difference in the 
boundary conditions which apply at the defect-superconductor (D-S) boundary for inclusion-like or roughness defects. 
Indeed, our calculations show, that for defects of the first type (with p < 0) the equality ^ ~ holds at the D-S 
boundary; that is in fact equivalent to effectively employing the N-S boundary condition at the D-S boundary. Note 
that the latter feature holds specifically for defects with p < 0; in the opposite case (p > 0) the effective boundary 
condition should be formulated anew. On the contrary, for defects of the second type we employ the I-S boundary 
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conditions (i.e. those between an insulator and superconductor) at the D-S boundary. 

To study the role of boundary conditions (N-S -type versus I-S -type) in more detail we have performed a series 
of numerical calculations in order to find H en for a k = 2 superconductor with ideal surface. We found out that 
the value of H en — 0.36i? C 2 (for N-S case) is lower than that for the I-S case (H en = 0AIH C 2 ). Thus, the type of 
boundary conditions appreciably affects the value of H en . 

The increase of H en found in the case of the inclusion defect with < p(r) < 1 (in comparison with the reference 
case p(r) = 0) is explained as follows. The greater p(r) the less vortices penetrate into a defect with T* < T c . This 
eventually decreases the value of II near the corners (as compared to its value in the case p(r) = 0); therefore a 
somewhat larger magnetic field is required for the kinematic momentum II to reach its critical value II*: ' . 

In addition to the above calculations we examined the distribution of the current density near defects of arbitrary 
rectangular shape in bulk superconductors. Our study was performed by numerical solution of the London equation 
applied to a superconducting slab in the Meissner state. It turns out that the current density j COI inside the corners 
of defects (or, more precisely, at a distance A/20 from the corners, the spacing of our grid) exceeds the value joo on 
the same slab surface far from the defect. Besides, we found that for defects whose size exceeds A, the ratio jcor/joo 
is practically constant, see table IV. 

According to our calculation the ratio jcor/joo is not very sensitive to the shape of the defect. For the case of a 
rectangular defect, it turns out that with the decrease of the defect width (at a fixed depth) the ratio j C or/ j<x increases 
somewhat. The above result qualitatively explains why H en saturates with increasing defect size, as reported in tables 
II, III. 

Thus, the surface defects are capable to appreciably suppress the barrier for the flux entry and exit in a supercon- 
ductor. However, this suppression is not complete, since the first vortex penetration field H en is still greater than 
the first critical field H c \- For example, for n = 2 one finds H en ~ 1A2H C \, for k = 5, H en ~ 1.93if c i, and for 
hi = 10, H en ~ 2.54H c i (we selected defects of size £ x 2A). These results are obtained for type-II narrow and long 
defects which, as it was shown recently 
according to our calculations, most effective are wide and long solitary defects of length greater than A and of width 
much greater than £ (the latter condition excludes the proximity effect at the D-S boundary). 



17 lq], suppress the barrier in a most effective way. In case of type-I defects, 



V. PHENOMENOLOGICAL MODEL OF AN EDGE/SURFACE BARRIER 

One of the main results of the present study is the weak dependence of the quantities II° dg0 (i7) and j° ds °(H) on the 
increasing magnetic field in a wide field range H en < H < H C 2, both with and without defects. This finding allows 
us to construct a phenomenological model of the surface barrier in bulk superconductors. Namely, we shall suppose 
that the vortices enter the sample when the current density at the surface/edge reaches some threshold value j s ; this 
quantity will be the main parameter of our model. For the case of an ideal surface, we will assume that the vortices 
exit from the superconductor when the current density on the surface vanishes. In the opposite case of a defective 
surface or edge, we assume that the vortices exit the superconductor through defects, thus leaving a vortex-free band 
(VFB) near the edges. Meissner currents flowing through those bands produce remanent magnetization during vortex 
exit. The width of the VFB, which depends on the surface quality, represents another parameter of the model. 

We consider this model taking as an example a wide slab free of surface defects. The calculations will be made in 
the framework of the nonlocal model of the critical state Although in this case we do not have bulk pinning, 

the term "critical state" obviously can be applied also in this case, since here there is a source of irreversibility - the 
surface barrier, which hinders not only the entry of vortices, but also their exit, as does bulk pinning. 

According to |19| ], the local induction B(y) satisfies the equation 

B-X 2 ^=n(y)^ 0l (19) 
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with boundary conditions B\ y= ± w /2 = H; here n(y) is the vortex density. In an increasing field the condition is added 
that the maximum current density (its absolute value) on the slab surface can not exceed j s . Since (19) contains two 
unknown functions B{y) and n(y), it is necessary to provide an additional condition, namely, the current density is 
zero where the density of vortices is finite and vice versa [gQ| . 

It is easy to show that at the initial stage of the increase of the magnetic field from zero up to H s = 
47rj s A/ctanh(w/2A) the density of vortices is equal to zero, i.e. the ordinary Meissner state takes place, and the 
solution of Eq. (19) reads 

fl(j/) = ffcosh(y/A)/cosh(u), (20) 

here v — w/2 A is a dimensionless parameter. With further increase of the field H the vortices start to penetrate into 
the slab, occupying the central area of the slab | y |< 6, where in the case v 3> 1, b is given by 

b^w/2 (l-±ln H + H * 



2v if - H s 

It is seen from the above equation that the vortices will occupy almost the entire sample [b(H ~ if, + 0) ~ w/2] 
practically at once, i.e. when if only slightly exceeds H s . This effect is a consequence of the short-range repulsion 
between Abrikosov vortices. 

The distribution of the magnetic field in the slab at given magnetic field is 



B(y) = 12/1 < 6 < (21) 

\ ff cosh^/2 -\y | )/\]-H. smh[(w/2-\y | )/A], b < \y\ < w/2 . 

Now we increase the magnetic field further up to some value H = Hq > H s and then decrease it down to zero. In this 
case, as follows from the model supplemented by Eq. (19), the central area occupied by the vortices expands until it 
reaches the surfaces of a slab; this happens at the exit field H = ff CX i t = (2b(Ho)/w) y/ if q — iff ~ a/ Hq — iff . In 
the range of fields if cx it < H < Hq the distribution of the magnetic field in a slab is as follows, 

[if/ cosh[K2 - &0/A], | V | < h, 

\ if cosh[(6! - |y|)/A]/cosh[(6i -w/2)/\], h < \y\ <w/2, 

where b\{H) is determined by the equation 

2b 1 H/H CKit = wcosh[(w/2 - b x )/\] . 

Knowing B(y) at each value of the magnetic field H, one can calculate the magnetization curve from Eqs. (20)-(22), 

f - H/4tt, 0<H<H s ,H T, 

M(H) = \ (V&=H* - B)/4*, H s <H<H Q ,H t (23) 
I (Ifexit - H)/4tt, if exit <H<H ,Hl, 

I 0, <H< if exit , iit- 

Note that the above expressions for the vortex-dome boundaries b(H), b\(H) slightly differ from the analogous 
expressions obtained by |po| . However, this does not produce a significant difference in the magnetization M(H), 
since one has b ~ w/2, b\ ~ w/2 within the actual range of the magnetic field in the limit d> 1. 

In Fig. 1 the dependences M(H) obtained from the phenomenological model for a pin-free superconductor (expres- 
sion (23) - curve lb) and for that with defects (curve 2b) are given. It is seen that in the absence of defects, the 
dependence M(H), calculated within the continuum approach, and that calculated from the GLE differ only slightly. 
This coincidence is quite surprising since the applicability area of the continuum approach for a slab of width 25£ 
seems rather narrow. Moreover, the continuum approach is valid in the case k> 1 because in the nonlocal theory the 
averaging is over distances exceeding the distance between vortices, but smaller than A. The paramagnetic response 
(with positive magnetization) resulting from the Ginzburg-Landau equations apparently cannot be reproduced within 
the framework of the phenomenological model of an edge/surface barrier due to the above mentioned limitation. 
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VI. CONCLUSION AND DISCUSSION 



In this paper using both numerical and analytical solutions of the time-dependent Ginzburg-Landau equations, 
the vortex entry condition for bulk and thin-film type-II superconductors is studied. A universal "supervelocity 
criterion" for the vortex entry into type-II superconductors is formulated, which is shown to be more precise than the 
conventional "supercurrent criterion" . The role of surface defects on the vortex entry condition and also on the shape of 
the magnetization curve is established. Vortices start to nucleate and penetrate into superconductors through surface 
defects. The exit of vortices also occurs through defects. As a result, the magnetization curves of superconductors 
with an ideal surface and with an imperfect surface differ not only quantitatively but also qualitatively. On the basis 
of the obtained results a simple phcnomcnological model of an edge/surface barrier is suggested. The comparison of 
magnetization curves obtained from our analytical model with that obtained numerically from the GLE, demonstrates 
quite good qualitative and in some respect quantitive agreement, despite the simplicity of our phenomenological model. 

The obtained results for the vortex entry conditions into thin-film and bulk isotropic superconductors suggest a 
possible condition for the first vortex entry into layered and granular superconductors. Namely, in the case of layered 
superconductors (or superconducting multilayers) the currents flow only within the superconducting layers. Therefore, 
the vortex entry criterion for the supervelocity II should be reached in one layer. In the case of strongly anisotropic 
superconductors of the BiSCCO-type the thickness of a superconducting layer is d s — 3 A and the distance between 
layers is di = 12 A. It can be shown |^lj that the value of II at the edge will be larger than the analogous value 
for isotropic superconductors by a factor of U (d s + di)/d s ~ 2.3 for the same applied magnetic field and geometry 
of the sample. This factor reflects the renormalization of the London penetration depth in layered superconductors, 
A' = \y/ {d s + di)/d s pl| , p2| . As a result, the first vortex entry field will be reduced by a factor yj {d s + dijjd~ s as 
compared to isotropic superconductors. However, since the quantity A' (not A) can be directly measured in layered 
superconductors, the effect of the layered structure is automatically included in the recorded value of H en . In the 
more topical case of artificial superconducting multilayers p2| , the value of A inside the layer is well determined. This 
allows us to estimate both A' and H en in such systems, see e.g. |21||). 

In the case of granular superconductors we have a superconducting medium containing nonsuperconducting bridges, 
i.e. type-I defects in our nomenclature. As can be seen from Fig. 7, the entry field of a J-vortex into a defect can be 
much lower than the field H c \. If the defect on the surface is connected to defects deep inside the superconductor, 
a J-vortex can move deep into the superconductor through defects. In this case the defect may play a double role: 
both as a pinning center and as a channel along which the motion of a J-vortex occurs. Taking into account these 
properties of granular superconductors we conclude that the process of vortex entry into granular superconductors 
differs qualitatively from above considered cases. Since a GL treatment, being too complex, is not well suited to 
describe the magnetic properties of Josephson-like granular systems, some simplified approach should be chosen, for 
example similar to that used by |p3| , f24f . 
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TABLE I. Dependence of the nucleation time t en and of the number of entering vortices N (at one side of the slab) on the 
applied magnetic field H . 



H/H c2 


0.4250 


0.4170 


0.4110 


0.4060 


0.4050 


0.4047 


nj dgc 


1.071 


1.051 


1.036 


1.023 


1.021 


1.020 


N 


8 


7 


6 


5 


5 


5 


ten / T 


30 


35 


55 


120 


280 


780 
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TABLE II. Vortex entry field H en for various inclusion defects . 



Defect size 




2£ x 2£ 


3£x3£ 


4£ x 4£ 


6£x6C 


8^x8? 


hI„/h c2 


0.38 


0.30 


0.27 


0.28 


0.28 


0.28 



TABLE III. Dependence of the vortex entry field H en on the size of geometrical defects. 



Defect size 


2£ x 2£ 


3£ x 3£ 


4£ x 4£ 


5£ x 5£ 


6£ x 6£ 


7£ x 7£ 


11? x 11^ 


hH/h c2 


0.35 


0.35 


0.34 


0.34 


0.34 


0.33 


0.33 



TABLE IV. Dependence of j CO i/joo on the size of defects of square and rectangular shape. 



Defect size 


A/4 x A/4 


A/4 x A/2 


A/4 x A 


A/4 x 1.5A 


A/4 x 2A 


A/4 x 2.5A 


jcor /joo 


1.66 


2.03 


2.37 


2.53 


2.61 


2.66 


Defect size 


A/4 x A/4 


A/2 x A/2 


A x A 


1.5A x 1.5A 


2A x 2A 


2.5A x 2.5A 


j cor /joo 


1.66 


1.88 


2.08 


2.18 


2.23 


2.25 



Figure captions. 

Fig. 1. Magnetization curve of a bulk superconductor in the absence (curves la, lb) and in the presence of surface 
defects (curves 2a, 2b). Curves la and 2a are obtained by numerical solution of the TDGLE, and curves lb and 2b 
from a phenomenological model of the surface barrier (see Section V). 

Fig. 2. Time evolution of the vortex nucleation process at t = 15t (a), t = 24t (b), t = 30r (c), and t = 180t (d) 
after turning on a magnetic field H = 0A2H c2 . 

Fig. 3. The dependence of the length-average of n° dgc on the applied magnetic field. Solid curve: superconductor 
without defects. Dashed curve: superconductor with surface defects of size 2£ x 2£. 

Fig. 4. The dependence of the length-average of j^ dge on the applied magnetic field. Solid curve: superconductor 
without defects. Dashed curve: superconductor with surface defects of size 2£ x 2£. 

Fig. 5. The dependence of II cr = P 1 / 2 on the wave-number k of the disturbance q x obtained from Eq. (10) for 
different values of magnetic field H and GL parameter k: (1) n = 2, (2) k = 5, (3) n = 10, (4) n = 100. (a) H = H s , 
(b) H = H c2 . 

Fig. 6. Contour lines of the order parameter in a bulk superconductor with type-I surface defects of size 4£ x 4£. 
Magnetic field H = O.QH c2 . The parameters of the superconductor are given in the text. 

Fig. 7. The dependence of the length-average of n° dgc on the applied magnetic field < H < H en in the presence 
of surface defects. Curves 1-5 correspond to defects of size 7£ x 7£, 6£ x 6£, 5£ x 5£, 4£ x 4£, and 3£ x 3£, respectively. 
Curve 6 corresponds to a superconductor without defects. 
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Fig.l D.Yu. Vodolazov et.al. 




Fig. 2 D.Yu. Vodolazov et.al. 




Fig. 3 D.Yu. Vodolazov et.al. 




Fig. 4 D.Yu. Vodolazov et.al. 




Fig. 5 D.Yu. Vodolazov et.al. 



. 6 D.Yu. Vodolazov et.al. 




Fig. 7 D.Yu. Vodolazov et.al. 



